Binder 0#
Bienvenue dans ce chapitre interactif utilisant ThebeLab !
Exemple de code interactif#
import matplotlib.pyplot as plt
import numpy as np
def plot_sine(frequency=1):
x = np.linspace(0, 30, 1000)
y = np.sin(frequency * x)
plt.figure(figsize=(8, 4))
plt.plot(x, y)
plt.title(f'Sine Wave with Frequency {frequency}')
plt.xlabel('x')
plt.ylabel('sin(x)')
plt.show()
plot_sine()
try:
import dolfinx
print("dolfinx importé avec succès")
except ImportError as e:
print(f"Erreur lors de l'importation de dolfinx : {e}")
import sys
print(f"Python path : {sys.path}")
dolfinx importé avec succès
Problème torsion d’arbres#
Librairies#
# Importation des librairies
import dolfinx # Module principal de FEniCSx pour le calcul par éléments finis
from dolfinx import mesh, fem, plot, io, default_scalar_type # Sous-modules FEniCSx essentiels
from dolfinx.fem.petsc import LinearProblem # Résolution de systèmes linéaires
from mpi4py import MPI # Interface pour le calcul parallèle
import ufl # Unified Form Language pour les formulations variationnelles
import numpy as np # Calcul numérique efficace
import matplotlib.pyplot as plt
import math # Module de fonctions mathématiques standard de Python
import pyvista as pv # Visualisation 3D scientifique
import gmsh # Générateur de maillage 3D
import meshio # Lecture/écriture de différents formats de maillage
from dolfinx.io import XDMFFile # Gestion des fichiers XDMF pour données volumineuses
import dolfinx.fem as fem # Module FEniCSx pour la méthode des éléments finis
import dolfinx.mesh as mesh # Module FEniCSx pour la gestion avancée des maillages
from petsc4py import PETSc # Suite de solveurs numériques pour problèmes à grande échelle
import panel as pn
from myst_nb import glue
import ipywidgets as widgets
from IPython.display import display
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
Géométrie#
Définition de la géométrie et du maillage#
# Initialisation du plotteur PyVista en dehors de la fonction
p = pv.Plotter(notebook=True)
p.set_background("grey")
# Fonction pour créer et visualiser la géométrie
def create_ellipse(rho_x, rho_y, h=1.0, lc=0.02):
global domain
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
gmsh.model.add("ellipse")
# Création de l'ellipse
ellipse = gmsh.model.occ.addEllipse(0, 0, 0, rho_x, rho_y)
curve_loop = gmsh.model.occ.addCurveLoop([ellipse])
surface = gmsh.model.occ.addPlaneSurface([curve_loop])
# Extrusion pour obtenir un cylindre elliptique
gmsh.model.occ.extrude([(2, surface)], 0, 0, h)
gmsh.model.occ.synchronize()
# Configuration du maillage
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", lc)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)
gmsh.model.mesh.generate(3)
# Sauvegarde et conversion du maillage
gmsh.write("ellipse.msh")
msh = meshio.read("ellipse.msh")
meshio.write("ellipse.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells_dict.get("tetra", [])}))
gmsh.finalize()
# Lecture du maillage avec FEniCSx
with XDMFFile(MPI.COMM_WORLD, "ellipse.xdmf", "r") as xdmf_file:
domain = xdmf_file.read_mesh(name="Grid")
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
# Conversion du maillage pour la visualisation avec PyVista
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(domain)
u_grid = pv.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
# Effacer la scène précédente et ajouter le nouveau maillage
p.clear()
# Ajout du maillage à la scène de visualisation
p.add_mesh(u_grid,
show_edges=True, # Affiche les arêtes du maillage
scalar_bar_args={
"title": "u", # Titre de la barre de couleur
"title_font_size": 24,
"label_font_size": 22,
"shadow": True,
"italic": True,
"font_family": "arial",
"vertical": False # Orientation horizontale de la barre de couleur
})
# Ajout d'un titre à la visualisation
p.add_text("Cylindre Mesh", font_size=24, color="black", position="upper_edge")
# Ajout des limites de la boîte englobante
p.show_bounds(color="black")
# Ajout des axes de coordonnées
p.add_axes(color="black")
# Définition de la couleur de fond
p.set_background("grey")
# Mise à jour de la scène
p.show()
return rho_x, rho_y, h
# Widgets interactifs pour a, b, et h
rho_x_slider = widgets.FloatSlider(value=0.2, min=0.1, max=1.0, step=0.01, description='Demi-grand axe')
rho_y_slider = widgets.FloatSlider(value=0.2, min=0.1, max=1.0, step=0.01, description='Demi-petit axe')
h_slider = widgets.FloatSlider(value=1.0, min=0.5, max=3.0, step=0.1, description='Hauteur')
# Interface utilisateur et fonction interactive
ui = widgets.VBox([rho_x_slider, rho_y_slider, h_slider])
out = widgets.interactive_output(create_ellipse, {'rho_x': rho_x_slider, 'rho_y': rho_y_slider, 'h': h_slider})
display(ui, out)
# Fonction pour récupérer les valeurs actuelles des curseurs
def get_slider_values():
current_rho_x = rho_x_slider.value
current_rho_y = rho_y_slider.value
current_h = h_slider.value
return current_rho_x, current_rho_y, current_h
# Exemple d'utilisation : récupération des valeurs
rho_x, rho_y, h = get_slider_values()
rho = np.array([rho_x, rho_y]) # Demi-axes dans les directions x et y
Unhandled exception
Traceback (most recent call last):
File "/opt/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/aiohttp/web_protocol.py", line 531, in start
resp, reset = await task
^^^^^^^^^^
File "/opt/anaconda3/envs/fenicsx-env/lib/python3.12/asyncio/futures.py", line 287, in __await__
yield self # This tells Task to wait for completion.
^^^^^^^^^^
File "/opt/anaconda3/envs/fenicsx-env/lib/python3.12/asyncio/tasks.py", line 385, in __wakeup
future.result()
File "/opt/anaconda3/envs/fenicsx-env/lib/python3.12/asyncio/futures.py", line 203, in result
raise self._exception.with_traceback(self._exception_tb)
File "/opt/anaconda3/envs/fenicsx-env/lib/python3.12/asyncio/tasks.py", line 314, in __step_run_and_handle_result
result = coro.send(None)
^^^^^^^^^^^^^^^
File "/opt/anaconda3/envs/fenicsx-env/lib/python3.12/site-packages/aiohttp/web_protocol.py", line 448, in _handle_request
assert self._request_handler is not None
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AssertionError
interactive_panel = pn.pane.VTK(p.ren_win, width=755, height=400)
from myst_nb import glue
glue("img", interactive_panel)
Configuration du problème#
Définition de l’espace fonctionnel#
# Définition de l'espace fonctionnel pour le problème
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
Définition des frontière du domaine#
# Définition des fonctions pour identifier les différentes faces du cylindre
def down_face(x):
"""Identifie la face inférieure du cylindre (z = 0)"""
return np.isclose(x[2], 0)
def top_face(x):
"""Identifie la face supérieure du cylindre (z = h)"""
return np.isclose(x[2], h)
def lateral_face(x):
"""
Identifie la surface latérale du cylindre elliptique
Utilise l'équation de l'ellipse : x^2/a^2 + y^2/b^2 = 1
"""
tolerance = 1e-5 # Tolérance pour la comparaison numérique
return (np.isclose((x[0]**2 / rho[0]**2 + x[1]**2 / rho[1]**2), 1.0, atol=tolerance)) & (0 <= x[2]) & (x[2] <= h)
# Dimension des facettes (2D pour un domaine 3D)
fdim = domain.topology.dim - 1
# Localisation des entités (facettes) correspondant à chaque face
Sigma_l = mesh.locate_entities(domain, fdim, lateral_face)
Sigma_0 = mesh.locate_entities(domain, fdim, down_face)
Sigma_h = mesh.locate_entities(domain, fdim, top_face)
# Préparation des marqueurs de facettes pour les conditions aux limites
marked_facets = np.hstack([Sigma_l, Sigma_0, Sigma_h])
marked_values = np.hstack([
np.full_like(Sigma_l, 1), # Marqueur 1 pour la face latérale
np.full_like(Sigma_0, 2), # Marqueur 2 pour la face inférieure
np.full_like(Sigma_h, 3) # Marqueur 3 pour la face supérieure
])
# Tri et création des tags de maillage
sorted_indices = np.argsort(marked_facets)
facet_tag = dolfinx.mesh.meshtags(domain, fdim,
marked_facets[sorted_indices],
marked_values[sorted_indices])
# Définition de la mesure pour les facettes marquées
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag)
# Création d'un dictionnaire pour mapper les noms aux valeurs numériques
surface_ids = {'Sl': 1, 'S0': 2, 'Sh': 3}
Visualisation des frontière du domaine#
# Définition d'une fonction pour appliquer des marqueurs aux facettes
def apply_marker(boundary_facets, marker_array, marker_value):
"""
Applique un marqueur spécifique aux facettes de la frontière.
Args:
boundary_facets: Indices des facettes de la frontière
marker_array: Tableau des marqueurs à mettre à jour
marker_value: Valeur du marqueur à appliquer
"""
# Filtrer les facettes pour ne garder que celles appartenant au domaine local
boundary_facets = boundary_facets[boundary_facets < num_cells_local]
# Appliquer le marqueur
marker_array[boundary_facets] = marker_value
# Obtenir le nombre de cellules locales (important pour le calcul parallèle)
num_cells_local = domain.topology.index_map(fdim).size_local
# Initialiser les tableaux de marqueurs pour chaque type de frontière
# On crée 4 tableaux : 3 pour chaque face et 1 pour la combinaison
markers = [np.zeros(num_cells_local, dtype=np.int32) for _ in range(4)]
# Appliquer les marqueurs spécifiques à chaque face
apply_marker(Sigma_0, markers[0], 1) # Face inférieure (z = 0)
apply_marker(Sigma_h, markers[1], 2) # Face supérieure (z = h)
apply_marker(Sigma_l, markers[2], 3) # Surface latérale
# Calculer le marqueur combiné
markers[3] = markers[0] + markers[1] + markers[2]
# Ajouter un marqueur spécifique (4) pour les cellules non marquées
# Cela peut être utile pour identifier le volume intérieur, par exemple
markers[3][markers[3] == 0] = 4
# Créer la connectivité topologique pour les facettes
# Ceci est nécessaire pour certaines opérations sur le maillage
domain.topology.create_connectivity(fdim, fdim)
# Obtenir les données de maillage au format compatible avec PyVista
topology, cell_types, x = dolfinx.plot.vtk_mesh(domain, fdim, np.arange(num_cells_local, dtype=np.int32))
# Création du maillage PyVista à partir des données FEniCSx
grid = pv.UnstructuredGrid(topology, cell_types, x)
def add_plot(ax, marker, color, title, threshold_min):
"""
Fonction pour ajouter des maillages à une sous-fenêtre de visualisation.
Args:
ax: Axe de la sous-fenêtre PyVista
marker: Tableau des marqueurs pour les cellules
color: Couleur pour les cellules filtrées
title: Titre de la sous-fenêtre
threshold_min: Valeur minimale pour le filtrage des cellules
"""
# Mise à jour des données de cellule du maillage avec les marqueurs
grid.cell_data["Marker"] = marker
grid.set_active_scalars("Marker")
# Ajout du maillage complet avec les marqueurs
ax.add_mesh(grid, show_edges=True, color="cyan", scalar_bar_args={
"title": "Boundary Marker",
"title_font_size": 24,
"label_font_size": 22,
"shadow": True,
"italic": True,
"font_family": "arial",
"vertical": False
})
# Application d'un filtre basé sur le seuil pour isoler certaines cellules
grid_filter = grid.threshold(threshold_min, scalars='Marker')
ax.add_mesh(grid_filter, color=color, show_edges=True)
# Ajout du titre et des axes à la sous-fenêtre
ax.add_text(title, font_size=12, color="black", position="upper_edge")
ax.add_axes(color="black")
# Configuration de la visualisation PyVista avec 4 sous-fenêtres (2x2)
pl = pv.Plotter(shape=(2, 2))
# Sous-fenêtre 1 : Affichage de toutes les frontières
pl.subplot(0, 0)
add_plot(pl, markers[3], "red", "All Boundaries", threshold_min=0.5)
# Sous-fenêtre 2 : Affichage de la frontière inférieure
pl.subplot(0, 1)
add_plot(pl, markers[0], "red", "Down Boundary", threshold_min=0.5)
# Sous-fenêtre 3 : Affichage de la frontière supérieure
pl.subplot(1, 0)
add_plot(pl, markers[1], "red", "Top Boundary", threshold_min=1.5)
# Sous-fenêtre 4 : Affichage de la frontière latérale
pl.subplot(1, 1)
add_plot(pl, markers[2], "red", "Lateral Boundary", threshold_min=2.5)
# Configuration finale et affichage
pl.set_background("grey")
pl.show()
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(pl.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)
Conditions aux limites#
# Définition des conditions aux limites de Dirichlet (déplacement nul)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)
# Application des conditions aux limites sur chaque face
bc0_S0 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_0), V) # Face inférieure
bc0_Sh = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_h), V) # Face supérieure
bc0_Sl = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_l), V) # Surface latérale
Propriétés matériau#
import ipywidgets as widgets
from IPython.display import display, clear_output
from ipywidgets import HBox, VBox
# Widgets pour les constantes élastiques en GPa
mu_input = widgets.FloatText(value=80, description="μ = ", step=1)
mu_label = widgets.Label(value="GPa")
lambda_input = widgets.FloatText(value=120, description="λ = ", step=1)
lambda_label = widgets.Label(value="GPa")
# Organisation avec les labels d'unité
mu_box = HBox([mu_input, mu_label])
lambda_box = HBox([lambda_input, lambda_label])
# Fonction pour obtenir les constantes élastiques en Pa
def display_elastic_constants(change):
mu = mu_input.value * 1e9 # Conversion en Pa
lambda_ = lambda_input.value * 1e9 # Conversion en Pa
# Afficher les valeurs en Pa
clear_output(wait=True) # Clear previous output
return mu, lambda_ # Retourner les valeurs en Pa
# Liaison de la fonction d'affichage aux événements de changement de valeur des widgets
mu_input.observe(display_elastic_constants, names='value')
lambda_input.observe(display_elastic_constants, names='value')
# Affichage des widgets
ui_elastic = VBox([mu_box, lambda_box])
display(ui_elastic)
mu, lambda_ = display_elastic_constants('value')
print(f"Module de cisaillement (μ) en Pa: {mu:.2f}")
print(f"Premier paramètre de Lamé (λ) en Pa: {lambda_:.2f}")
Module de cisaillement (μ) en Pa: 80000000000.00
Premier paramètre de Lamé (λ) en Pa: 120000000000.00
Définition des chargements#
import ipywidgets as widgets
from IPython.display import display
import numpy as np
def create_angle_widget():
# Création du widget curseur
angle_slider = widgets.FloatSlider(
value=2.0, # Valeur initiale à 2 degrés
min=0,
max=360,
step=1,
description='Angle (°):',
disabled=False,
continuous_update=True,
orientation='horizontal',
readout=True,
readout_format='.0f',
)
# Widget de sortie pour afficher les valeurs
output = widgets.Output()
# Fonction pour mettre à jour l'affichage
def update_angle(change):
ang_degrees = change['new']
ang_radians = np.radians(ang_degrees)
with output:
output.clear_output(wait=True)
print(f"Angle : {ang_degrees:.0f}° ({ang_radians:.2f} radians)")
return ang_radians
# Attacher la fonction de mise à jour au widget
angle_slider.observe(update_angle, names='value')
# Afficher le widget et la sortie
display(widgets.VBox([angle_slider, output]))
# Initialiser l'affichage
update_angle({'new': angle_slider.value})
return angle_slider
# Créer et afficher le widget
angle_widget = create_angle_widget()
# Fonction pour obtenir la valeur actuelle en radians
def get_angle_radians():
return np.radians(angle_widget.value)
# Exemple d'utilisation
# Obtention de l'angle de torsion en radians à partir du widget
ang_radians = get_angle_radians()
# Calcul de alpha (taux de torsion)
# alpha représente le taux de rotation par unité de longueur
# Il est calculé en divisant l'angle de torsion par la longueur du cylindre (h)
# Unités : radians / mètre (rad/m)
alpha = ang_radians / h
print(f"Taux de torsion (α) = {alpha:.4f} rad/m")
# Calcul du rayon équivalent du cylindre elliptique
# On utilise la moyenne quadratique des demi-axes pour obtenir un rayon équivalent
# Cette formule donne une bonne approximation pour les calculs de torsion
R = np.sqrt((rho[0]**2 + rho[1]**2) / 2)
print(f"Rayon équivalent (R) = {R:.4f} m")
# Calcul de la rigidité en torsion C
pi = math.pi # Nombre Pi
# La formule pour C est basée sur la théorie de la torsion pour un cylindre circulaire
# mu : module de cisaillement du matériau
# R^4 : le rayon à la puissance 4 représente l'influence de la géométrie
# alpha : taux de torsion calculé précédemment
# Le facteur 1/2 est spécifique à la formule de rigidité en torsion
C = mu * pi * (R ** 4) / 2 * alpha
print(f"Rigidité en torsion (C) = {C:.4e} N·m²")
Taux de torsion (α) = 0.0349 rad/m
Rayon équivalent (R) = 0.2000 m
Rigidité en torsion (C) = 7.0184e+06 N·m²
# Calcul de la constante A pour le champ de torsion
A = 2 * C / (pi * R**4) # A est un scalaire constant positif
print(f"Densité de contrainte de torsion (A_τ) = {A:.4e} N/m³")
# Définition des coordonnées spatiales
x = ufl.SpatialCoordinate(domain)
# Définition des composantes du champ de torsion (vecteur des Contraintes de Cauchy)
T1 = -A * x[1] # Composante x du champ de torsion / Contrainte de cisaillement dans la direction x
T2 = A * x[0] # Composante y du champ de torsion / Contrainte de cisaillement dans la direction y
T3 = ufl.as_ufl(0.0) # Composante z nulle / Pas de contrainte dans la direction z (axe de torsion)
# Combinaison en un champ vectoriel de torsion
T = ufl.as_vector([T1, T2, T3])
print("Vecteur de contrainte de Cauchy défini pour la torsion (N/m²)")
Densité de contrainte de torsion (A_τ) = 2.7925e+09 N/m³
Vecteur de contrainte de Cauchy défini pour la torsion (N/m²)
# Définition des tenseurs de déformation et de contrainte
def epsilon(u):
"""Tenseur de déformation linéarisé"""
return ufl.sym(ufl.grad(u))
def sigma(u):
"""Tenseur des contraintes (loi de Hooke)"""
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)
# Définition des fonctions d'essai et de test
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# Définition de la force volumique (nulle dans ce cas)
f = fem.Constant(domain, default_scalar_type((0, 0, 0)))
# Formulation variationnelle du problème
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx # Forme bilinéaire
# Créez un widget Dropdown
surface_selector = widgets.Dropdown(
options=[('Surface latérale Sl', 'Sl'), ('Surface inférieure S0', 'S0'), ('Surface supérieure Sh', 'Sh')],
value='Sh',
description='Surface:',
disabled=False,
)
# Créez un widget de sortie pour afficher les messages
output = widgets.Output()
# Fonction pour mettre à jour la forme linéaire
def update_linear_form(change):
global L
selected_surface = change['new']
# Mise à jour de la forme linéaire
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds(surface_ids[selected_surface])
# Mise à jour de l'affichage
with output:
output.clear_output(wait=True)
print(f"Torsion appliquée sur {selected_surface}")
# Attacher la fonction de mise à jour au widget
surface_selector.observe(update_linear_form, names='value')
# Afficher le widget et la sortie
display(widgets.VBox([surface_selector, output]))
# Initialiser l'affichage
update_linear_form({'new': surface_selector.value})
# Créez un widget Dropdown pour la sélection de la condition aux limites
bc_selector = widgets.Dropdown(
options=[
('Aucune', 'none'),
('Face inférieure S0', 'S0'),
('Face supérieure Sh', 'Sh'),
('Surface latérale Sl', 'Sl')
],
value='S0', # Valeur par défaut
description='Encastrement:',
disabled=False,
)
# Dictionnaire pour mapper les sélections aux conditions aux limites
bc_dict = {
'none': [],
'S0': [bc0_S0],
'Sh': [bc0_Sh],
'Sl': [bc0_Sl]
}
# Créez un widget de sortie pour afficher les messages
output = widgets.Output()
# Fonction pour mettre à jour le problème
def update_problem(change):
global problem
selected_bc = change['new']
problem = LinearProblem(a, L, bcs=bc_dict[selected_bc],
petsc_options={"ksp_type": "cg", "pc_type": "jacobi"})
# Mise à jour de l'affichage
with output:
output.clear_output(wait=True)
print(f"Encastrement appliquée sur {selected_bc}")
# Attacher la fonction de mise à jour au widget
bc_selector.observe(update_problem, names='value')
# Afficher le widget
display(bc_selector, output)
# Initialisation du problème avec la valeur par défaut
update_problem({'new': bc_selector.value})
uh = problem.solve()
# Création du maillage pour PyVista basé sur les coordonnées des dofs
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(domain)
# Créez la grille PyVista et ajoutez les valeurs des dofs à la grille
u_grid = pv.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
# Attach vector values to grid and warp grid by vector
u_grid["u"] = uh.x.array.reshape((u_geometry.shape[0], 3))
# Visualisation
p = pv.Plotter()
p.add_mesh(u_grid, show_edges=False, scalar_bar_args={
"title": "u",
"title_font_size": 24,
"label_font_size": 22,
"shadow": False,
"italic": True,
"font_family": "arial",
"vertical": False
})
p.add_text("Déplacements", font_size=12, color="black", position="upper_edge")
p.add_axes(color="black")
p.set_background("grey")
p.show()
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)
Visualisation des rotations principales#
def omega(u):
"""Calcule le tenseur de rotation (partie antisymétrique du gradient de déplacement)"""
return 0.5 * (ufl.grad(u).T - ufl.grad(u))
# Calcul du tenseur de rotation
theta = omega(uh)
# Extraction des composantes du tenseur de rotation
theta_xy, theta_xz, theta_yz = theta[0, 1], theta[0, 2], theta[1, 2]
# Définition de l'espace fonctionnel pour les rotations (éléments discontinus d'ordre 0)
V_omega = fem.functionspace(domain, ("DG", 0))
def interpolate_rotation(rotation_component):
"""Interpole une composante de rotation sur l'espace fonctionnel"""
expr = fem.Expression(rotation_component, V_omega.element.interpolation_points())
rotation = fem.Function(V_omega)
rotation.interpolate(expr)
return rotation
# Interpolation des composantes de rotation
rotation_x = interpolate_rotation(theta_yz)
rotation_y = interpolate_rotation(-theta_xz) # Note le signe négatif
rotation_z = interpolate_rotation(theta_xy)
# Calcul de la norme de rotation
norm_theta = ufl.sqrt(theta_xy**2 + theta_xz**2 + theta_yz**2)
rotation_norm = interpolate_rotation(norm_theta)
# Configuration de la visualisation
pl = pv.Plotter(shape=(2, 2))
dargs = dict(cmap="jet", show_scalar_bar=False)
warped = u_grid.warp_by_vector("u", factor=0) # Pas de déformation pour visualiser les déformations
def add_subplot(row, col, data, title):
"""Ajoute un sous-plot pour une composante de rotation"""
pl.subplot(row, col)
warped.cell_data[title] = data.vector.array
warped.set_active_scalars(title)
pl.add_mesh(warped, **dargs)
pl.add_axes(color="k")
pl.add_text(title, color='k', font_size=12)
# Création des sous-plots
add_subplot(0, 0, rotation_norm, "Norme de rotation")
add_subplot(0, 1, rotation_x, "Rotation RX")
add_subplot(1, 0, rotation_y, "Rotation RY")
add_subplot(1, 1, rotation_z, "Rotation RZ")
# Configuration finale
pl.link_views()
pl.camera_position = 'iso'
pl.background_color = 'grey'
# Ajout d'une barre de couleur globale
pl.add_scalar_bar(title="Rotation (rad)", n_labels=5, label_font_size=10, title_font_size=12,
position_x=0.25, position_y=0.05, width=0.5)
# Affichage
pl.show()
Visualisation des déplacements amplifiés#
# Création d'un objet Plotter PyVista pour la visualisation 3D
p = pv.Plotter()
# Application de la déformation au maillage
# Le facteur 0.04 amplifie la déformation pour une meilleure visibilité
warped = u_grid.warp_by_vector("u", factor=0.07)
# Ajout du maillage déformé à la scène
p.add_mesh(warped,
show_edges=False, # Affiche les arêtes du maillage
scalar_bar_args={
"title": "Déplacement", # Titre de la barre de couleur
"title_font_size": 24,
"label_font_size": 22,
#"shadow": True,
"italic": True,
"font_family": "arial",
"vertical": False, # Orientation horizontale de la barre de couleur
"n_labels": 5, # Nombre d'étiquettes sur la barre de couleur
#"fmt": "%.2e" # Format scientifique pour les valeurs
})
# Ajout d'un titre à la visualisation
p.add_text("Déplacements amplifiés", font_size=12, color="k", position="upper_edge")
# Ajout des axes de coordonnées
p.add_axes(xlabel='X', ylabel='Y', zlabel='Z', color="k")
# Configuration de l'arrière-plan et de l'éclairage
p.set_background("grey")
#p.enable_shadows() # Ajoute des ombres pour améliorer la perception 3D
# Configuration de la caméra pour une vue optimale
p.camera_position = [(3, 3, 2), (0, 0, 0.5), (0, 0, 1)]
# Affichage de la visualisation
p.show()
Visualisation des déplacements normalisés et selon les axes principaux#
# Définition des arguments communs pour l'affichage des maillages
dargs = dict(
scalars="u", # Utilise le champ de déplacement 'u' pour la coloration
cmap="jet", # Utilise la palette de couleurs 'jet'
show_scalar_bar=False, # Ne pas afficher la barre de couleur pour chaque sous-plot
)
# Création d'un plotter PyVista avec une grille 2x2 de sous-plots
pl = pv.Plotter(shape=(2, 2))
# Sous-plot 0,0 : Déplacement normalisé (magnitude)
pl.subplot(0, 0)
pl.add_mesh(u_grid, **dargs)
pl.add_axes(color="black")
pl.add_text("Normalized Displacement", color='k', font_size=10)
# Sous-plot 0,1 : Déplacement selon X
pl.subplot(0, 1)
pl.add_mesh(u_grid.copy(), component=0, **dargs)
pl.add_axes(color="black")
pl.add_text("X Displacement", color='k', font_size=10)
# Sous-plot 1,0 : Déplacement selon Y
pl.subplot(1, 0)
pl.add_mesh(u_grid.copy(), component=1, **dargs)
pl.add_axes(color="black")
pl.add_text("Y Displacement", color='k', font_size=10)
# Sous-plot 1,1 : Déplacement selon Z
pl.subplot(1, 1)
pl.add_mesh(u_grid.copy(), component=2, **dargs)
pl.add_axes(color="black")
pl.add_text("Z Displacement", color='k', font_size=10)
# Configuration globale de la visualisation
pl.link_views() # Lie les vues des sous-plots pour une navigation synchronisée
pl.camera_position = 'iso' # Vue isométrique
pl.background_color = 'grey' # Couleur de fond grise
# Ajout d'une barre de couleur globale
pl.add_scalar_bar(title="Displacement (m)", n_labels=5, label_font_size=10, title_font_size=12, position_x=0.25, position_y=0.05, width=0.5)
# Affichage de la visualisation
pl.show()
Visualisation des composantes du tenseur des déformations#
# Calcul du tenseur de déformation
eps = epsilon(uh)
# Extraction des composantes du tenseur de déformation
eps_xx, eps_xy, eps_xz = eps[0, 0], eps[0, 1], eps[0, 2]
eps_yy, eps_yz = eps[1, 1], eps[1, 2]
eps_zz = eps[2, 2]
# Création de la grille PyVista
eps_topology, eps_cell_types, eps_geometry = dolfinx.plot.vtk_mesh(domain)
eps_grid = pv.UnstructuredGrid(eps_topology, eps_cell_types, eps_geometry)
# Définition de l'espace fonctionnel pour les déformations (éléments discontinus d'ordre 0)
V_eps = fem.functionspace(domain, ("DG", 0))
# Fonction pour interpoler les composantes de déformation
def interpolate_strain(eps_component):
expr = fem.Expression(eps_component, V_eps.element.interpolation_points())
strain = fem.Function(V_eps)
strain.interpolate(expr)
return strain
# Interpolation des composantes de déformation
strain_components = {
"XX": interpolate_strain(eps_xx),
"YY": interpolate_strain(eps_yy),
"ZZ": interpolate_strain(eps_zz),
"YZ": interpolate_strain(eps_yz),
"XZ": interpolate_strain(eps_xz),
"XY": interpolate_strain(eps_xy)
}
# Configuration de la visualisation
pl = pv.Plotter(shape=(2, 3))
warped = u_grid.warp_by_vector("u", factor=0) # Pas de déformation pour visualiser les déformations
dargs = dict(cmap="jet", show_scalar_bar=False)
# Fonction pour ajouter un sous-plot
def add_subplot(row, col, component):
pl.subplot(row, col)
warped.cell_data[f"Epsilon_{component}"] = strain_components[component].vector.array
warped.set_active_scalars(f"Epsilon_{component}")
pl.add_mesh(warped, **dargs)
pl.add_axes(color="black")
pl.add_text(f"Epsilon {component}", color='k', font_size=12)
# Création des sous-plots pour chaque composante
for i, component in enumerate(["XX", "YY", "ZZ", "YZ", "XZ", "XY"]):
add_subplot(i // 3, i % 3, component)
# Configuration finale
pl.link_views()
pl.camera_position = 'iso'
pl.background_color = 'grey'
# Ajout d'une barre de couleur globale
pl.add_scalar_bar(title="Strain", n_labels=5, label_font_size=10, title_font_size=12, position_x=0.25, position_y=0.05, width=0.5)
# Affichage
pl.show()
Visualisation des composantes du tenseur des contraintes#
# Calculer le tenseur des contraintes
sig = sigma(uh) # sigma est une fonction définie précédemment pour calculer les contraintes
# Extraire les composantes du tenseur des contraintes
sig_xx, sig_xy, sig_xz = sig[0, 0], sig[0, 1], sig[0, 2]
sig_yy, sig_yz = sig[1, 1], sig[1, 2]
sig_zz = sig[2, 2]
# Définir l'espace fonctionnel pour les contraintes (éléments discontinus d'ordre 0)
V_sig = fem.functionspace(domain, ("DG", 0))
# Fonction pour interpoler les composantes de contrainte
def interpolate_stress(sig_component):
expr = fem.Expression(sig_component, V_sig.element.interpolation_points())
stress = fem.Function(V_sig)
stress.interpolate(expr)
return stress
# Interpoler chaque composante de contrainte
stress_xx = interpolate_stress(sig_xx)
stress_yy = interpolate_stress(sig_yy)
stress_zz = interpolate_stress(sig_zz)
stress_yz = interpolate_stress(sig_yz)
stress_xz = interpolate_stress(sig_xz)
stress_xy = interpolate_stress(sig_xy)
# Configuration de la visualisation
pl = pv.Plotter(shape=(2, 3)) # Créer un plotter avec une grille 2x3
dargs = dict(
cmap="jet",
show_scalar_bar=False,
)
# Fonction pour ajouter un sous-plot
def add_subplot(row, col, stress_component, title):
pl.subplot(row, col)
warped.cell_data[title] = stress_component.vector.array
warped.set_active_scalars(title)
pl.add_mesh(warped, **dargs)
pl.add_axes(color="k")
pl.add_text(title, color='k', font_size=12)
# Ajouter chaque composante de contrainte comme un sous-plot
add_subplot(0, 0, stress_xx, "Sigma XX")
add_subplot(0, 1, stress_yy, "Sigma YY")
add_subplot(0, 2, stress_zz, "Sigma ZZ")
add_subplot(1, 0, stress_yz, "Sigma YZ")
add_subplot(1, 1, stress_xz, "Sigma XZ")
add_subplot(1, 2, stress_xy, "Sigma XY")
# Configuration finale de la visualisation
pl.link_views() # Lier les vues pour une navigation synchronisée
pl.camera_position = 'iso' # Vue isométrique
pl.background_color = 'grey' # Couleur de fond grise
# Ajout d'une barre de couleur globale
pl.add_scalar_bar(title="Contrainte (Pa)", n_labels=5, label_font_size=10, title_font_size=12, position_x=0.25, position_y=0.05, width=0.5)
# Afficher la visualisation
pl.show()
Calcul des contraintes de von Mises#
s = sigma(uh) - 1. / 3 * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))
von_Mises = ufl.sqrt(3. / 2 * ufl.inner(s, s))
V_von_mises = fem.functionspace(domain, ("DG", 0))
stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
stresses = fem.Function(V_von_mises)
stresses.interpolate(stress_expr)
Visualisation des contraintes de von Mises#
# Créez une nouvelle instance de Plotter
# Ajout des données de contrainte de von Mises au maillage déformé
warped.cell_data["VonMises"] = stresses.vector.array
# Définition de "VonMises" comme scalaire actif pour la visualisation
warped.set_active_scalars("VonMises")
# Création d'une nouvelle instance de Plotter PyVista
p = pv.Plotter()
# Ajout du maillage déformé avec la contrainte de von Mises
p.add_mesh(warped,
cmap="jet", # Utilisation de la carte de couleurs "jet"
show_edges=False, # Ne pas afficher les arêtes du maillage
scalar_bar_args={
"title": "von Mises (MPa)", # Titre de la barre de couleur
"title_font_size": 24, # Taille de la police du titre
"label_font_size": 22, # Taille de la police des étiquettes
"shadow": True, # Ajout d'une ombre à la barre de couleur
"italic": True, # Texte en italique
"font_family": "arial", # Police Arial
"vertical": False # Barre de couleur horizontale
})
# Ajout d'un titre à la visualisation
p.add_text("Contraintes de von Mises", font_size=12, color="black", position="upper_edge")
# Ajout des axes de coordonnées en noir
p.add_axes(color="black")
# Configuration de l'arrière-plan en gris
p.set_background("grey")
# Configuration de la caméra pour une vue optimale
p.camera_position = [(3, 3, 2), (0, 0, 0.5), (0, 0, 1)]
# Affichage de la visualisation
p.show()
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)
Analyses des résultats#
Calcul des valeurs maximales#
# Calcul du moment d'inertie polaire J pour une section circulaire
J = pi * (R ** 4) / 2
# Calcul de l'angle de torsion théorique
theta_theoretical = C * h / (mu * J)
# Extraction de la valeur maximale de la rotation RZ
max_rotation_z = np.max(np.abs(rotation_z.x.array))
# Extraction de la valeur maximale de la norme de rotation
max_rotation_norm = np.max(rotation_norm.x.array)
# Conversion en degrés pour une interprétation plus intuitive
max_rotation_z_degrees = np.degrees(max_rotation_z)
# Calcul de l'angle de torsion théorique
theoretical_twist = ang_radians * h # ang est l'angle de torsion par unité de longueur, h est la hauteur du cylindre
# Comparaison pour RZ
relative_error_Z = abs(max_rotation_z - theoretical_twist) / theoretical_twist * 100
# Comparaison pour norme R
relative_error_N = (max_rotation_norm - theta_theoretical) / theta_theoretical * 100
print(np.max(rotation_z.x.array))
0.03471230845573148
# Calcul du moment d'inertie polaire J pour une section circulaire
J = pi * (R ** 4) / 2
# Calcul de l'angle de torsion théorique
theta_theoretical = C * h / (mu * J)
# Affichage des résultats
print(f"Moment d'inertie polaire J : {J:.6e} m^4")
print(f"Angle de torsion théorique θ : {theta_theoretical:.6f} radians")
print(f"Angle de torsion théorique θ : {np.degrees(theta_theoretical):.2f} degrés")
# Comparaison avec la rotation maximale calculée numériquement
if 'max_rotation_z' in locals():
relative_error = abs(max_rotation_z - theta_theoretical) / theta_theoretical * 100
print(f"\nComparaison avec la simulation numérique en RZ:")
print(f"Rotation Z maximale simulée : {max_rotation_z:.6f} radians")
print(f"Rotation Z maximale simulée : {np.degrees(max_rotation_z):.2f} degrés")
print(f"Erreur relative : {relative_error_Z:.2f}%")
print(f"\nComparaison avec la simulation numérique en norme:")
print(f"Norme de rotation maximale : {max_rotation_norm:.6f} radians")
print(f"Norme de rotation maximale : {np.degrees(max_rotation_norm):.2f} degrés")
print(f"Erreur relative : {relative_error_N:.2f}%")
else:
print("\nAttention : La rotation maximale simulée n'a pas été calculée précédemment.")
# Analyse des autres composantes
max_rotation_x = np.max(np.abs(rotation_x.x.array))
max_rotation_y = np.max(np.abs(rotation_y.x.array))
print("\nValeurs maximales des composantes de rotation :")
print(f"RX max : {np.degrees(max_rotation_x):.2f} degrés")
print(f"RX max : {max_rotation_x:.6f} radians")
print(f"RY max : {np.degrees(max_rotation_y):.2f} degrés")
print(f"RY max : {max_rotation_y:.6f} radians")
Moment d'inertie polaire J : 2.513274e-03 m^4
Angle de torsion théorique θ : 0.034907 radians
Angle de torsion théorique θ : 2.00 degrés
Comparaison avec la simulation numérique en RZ:
Rotation Z maximale simulée : 0.034712 radians
Rotation Z maximale simulée : 1.99 degrés
Erreur relative : 0.56%
Comparaison avec la simulation numérique en norme:
Norme de rotation maximale : 0.034734 radians
Norme de rotation maximale : 1.99 degrés
Erreur relative : -0.49%
Valeurs maximales des composantes de rotation :
RX max : 0.22 degrés
RX max : 0.003798 radians
RY max : 0.21 degrés
RY max : 0.003721 radians
# Extraction des valeurs maximales de rotation pour chaque composante
max_rotation_x = np.max(np.abs(rotation_x.x.array))
max_rotation_y = np.max(np.abs(rotation_y.x.array))
max_rotation_z = np.max(np.abs(rotation_z.x.array))
# Affichage des résultats
print("Valeurs maximales des composantes de rotation :")
print(f"RX max : {max_rotation_x:.6f} radians ({np.degrees(max_rotation_x):.2f} degrés)")
print(f"RY max : {max_rotation_y:.6f} radians ({np.degrees(max_rotation_y):.2f} degrés)")
print(f"RZ max : {max_rotation_z:.6f} radians ({np.degrees(max_rotation_z):.2f} degrés)")
# Calcul et affichage de la norme de rotation maximale
max_rotation_norm = np.sqrt(max_rotation_x**2 + max_rotation_y**2 + max_rotation_z**2)
print(f"\nNorme de rotation maximale : {max_rotation_norm:.6f} radians ({np.degrees(max_rotation_norm):.2f} degrés)")
# Comparaison avec l'angle de torsion théorique (si disponible)
if 'theta_theoretical' in locals():
print(f"\nAngle de torsion théorique : {theta_theoretical:.6f} radians ({np.degrees(theta_theoretical):.2f} degrés)")
# Calcul de l'erreur relative
relative_error = (max_rotation_z - theta_theoretical) / theta_theoretical * 100
print(f"Erreur relative (RZ max vs théorie) : {relative_error:.2f}%")
# Analyse de la contribution de chaque composante
total_rotation = max_rotation_x + max_rotation_y + max_rotation_z
print("\nContribution relative de chaque composante à la rotation totale :")
print(f"RX : {max_rotation_x / total_rotation * 100:.2f}%")
print(f"RY : {max_rotation_y / total_rotation * 100:.2f}%")
print(f"RZ : {max_rotation_z / total_rotation * 100:.2f}%")
Valeurs maximales des composantes de rotation :
RX max : 0.003798 radians (0.22 degrés)
RY max : 0.003721 radians (0.21 degrés)
RZ max : 0.034712 radians (1.99 degrés)
Norme de rotation maximale : 0.035117 radians (2.01 degrés)
Angle de torsion théorique : 0.034907 radians (2.00 degrés)
Erreur relative (RZ max vs théorie) : -0.56%
Contribution relative de chaque composante à la rotation totale :
RX : 8.99%
RY : 8.81%
RZ : 82.20%
import numpy as np
# Valeur analytique de la rotation maximale (ici, 0.035 radians)
theta_analytique = 0.035
# Valeur numérique récupérée
theta_numerique = max_rotation_z
# Calcul de l'erreur absolue
erreur_absolue = np.abs(theta_analytique - theta_numerique)
# Calcul de l'erreur relative
erreur_relative = erreur_absolue / np.abs(theta_analytique) * 100
# Affichage des résultats
print(f"Rotation analytique maximale : {theta_analytique} radians")
print(f"Rotation numérique maximale : {theta_numerique} radians")
print(f"Erreur absolue : {erreur_absolue} radians")
print(f"Erreur relative : {erreur_relative} %")
Rotation analytique maximale : 0.035 radians
Rotation numérique maximale : 0.03471230845573148 radians
Erreur absolue : 0.0002876915442685257 radians
Erreur relative : 0.8219758407672162 %